flowchart A[Define DAG] --> B(24 hours previous rain) A --> C(48 hours previous rain) B --> D(Define adjustment sets for each predictor) C --> D D --> E(Run all models satisfying\nthe back door criterion) E --> F(Average posterior probabilities) F --> G(Combine models in a single graph) style A fill:#44015466 style B fill:#3E4A894D style C fill:#3E4A894D style D fill:#26828E4D style E fill:#6DCD594D style F fill:#FDE7254D style G fill:#31688E4D
Statistical analysis using causal inference
Agalychnis lemur
Statistical analysis for the paper
In review. Environmental drivers of calling activity in the endangered species Lemur Leaf frog
Purpose
- Determine the adjustment sets that allow to infer a causal effect of environmental variables on vocal activity
- Evaluate causal effect of environmental factors on vocal activity of A. lemur with bayesian regression models
Instructions
The only input data needed is the file ‘acoustic_and_climatic_data_by_hour.csv’ which has been shared as supplementary data
Statistical models can take some time to fit, depending on computing power
Many output files are created and re-read in subsequent chunks
Analysis flowchart
This flowchart depicts the data analysis steps described in this report:
Load packages
The code below installs all the necessary packages to run the analyses described in the report:
Code
if (!require("sketchy", character.only = TRUE))
install.packages("sketchy")
pkgs <-
c(
"remotes",
"viridis",
"brms",
"cowplot",
"posterior",
"readxl",
"HDInterval",
"kableExtra",
"knitr",
"ggdist",
"ggplot2",
"lunar",
"cowplot",
"maRce10/brmsish",
"warbleR",
"ohun",
"dagitty",
"ggdag",
"tidybayes",
"pbapply"
)
# install/ load packages
sketchy::load_packages(packages = pkgs)
options("digits" = 6, "digits.secs" = 5, knitr.table.format = "html")
# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE, tidy = TRUE)
# set working directory as project directory or one directory above,
opts_knit$set(root.dir = "..")Custom functions
Here we create some functions to format data and model outputs:
Code
adjustment_set_formulas <- function(dag, exposure, required_variable, outcome, effect = "total",
type = "minimal", formula_parts = c(outcome), latent = NULL, remove = NULL, plot = TRUE,
...) {
if (plot)
gg <- ggdag_adjustment_set(.tdy_dag = tidy_dagitty(dag), exposure = exposure,
outcome = outcome, ...) + theme_dag()
temp_set <- adjustmentSets(x = dag, exposure = exposure, outcome = outcome, effect = effect,
type = type)
form_set <- lapply(temp_set, function(x) {
if (!is.null(remove))
x <- x[!x %in% remove]
form <- paste(formula_parts[1], " ~ ", exposure, " + ", paste(x, collapse = " + "),
if (length(formula_parts) == 2)
formula_parts[2] else NULL)
return(form)
})
form_set <- form_set[!duplicated(form_set)]
if (!is.null(latent))
for (i in latent) form_set <- form_set[!grepl(paste0(" ", i, " "), form_set)]
# form_set <- sapply(form_set, as.formula)
names(form_set) <- seq_along(form_set)
# add formula as attribute
attributes(form_set)$exposure.formula <- paste(formula_parts[1], " ~ ", exposure,
if (length(formula_parts) == 2)
formula_parts[2] else NULL)
if (plot)
return(gg) else return(form_set)
}
# Define a function to remove special characters
remove_special_chars <- function(text) {
# Replace special characters with hyphen
cleaned_text <- gsub("[^a-zA-Z0-9]+", "-", text)
# Remove leading and trailing hyphens
cleaned_text <- gsub("^-+|-+$", "", cleaned_text)
return(cleaned_text)
}
pa <- function(...) brms::posterior_average(...)
# to get average posterior values from models with different formulas
averaged_model <- function(formulas, data, model_call, ndraws = 1000, save.files = TRUE,
path = ".", suffix = NULL, cores = 1, name = NULL) {
if (dir.exists(file.path(path, name))) {
cat("Directory already existed. Attempting to fit missing models\n")
cat("Fitting models (step 1 out of 2) ...")
} else dir.create(path = file.path(path, name))
cat("Fitting models (step 1 out of 2) ...")
fit_list <- pblapply_brmsish_int(X = formulas, cl = cores, function(y) {
# make file name without special characters
mod_name <- paste0(remove_special_chars(as.character(y)), ".RDS")
if (save.files & !file.exists(file.path(path, mod_name))) {
cat("Fitting", y, "\n")
mc <- gsub(pattern = "formula", replacement = as.character(y), x = model_call)
mc <- parse(text = mc)
fit <- eval(mc)
if (save.files)
saveRDS(fit, file = file.path(path, mod_name))
} else {
cat("Reading", y, "(already existed)\n")
fit <- readRDS(file.path(path, mod_name))
}
return(fit)
})
if (length(formulas) > 1) {
cat("Averaging models (step 2 out of 2) ...")
average_call <- parse(text = paste("pa(", paste(paste0("fit_list[[", seq_along(fit_list),
"]]"), collapse = ", "), ", ndraws = ", ndraws, ")"))
# Evaluate the expression to create the call object
average_eval <- eval(average_call)
# add formula as attribute
attributes(average_eval)$averaged_fit_formulas <- formulas
rds_name <- if (is.null(suffix))
file.path(path, paste0(name, ".RDS")) else file.path(path, paste0(suffix, "_", name, ".RDS"))
if (save.files)
saveRDS(average_eval, file = rds_name)
# return draws from average models
return(average_eval)
} else cat("No model averaging conducted as a single formula was supplied")
}
# function to convert effect sizes to percentage of change
to_change_percentage <- function(x) {
(exp(x) - 1) * 100
}
# prints the results of the averaged models
draw_extended_summary <- function(draws, name = NULL, highlight = TRUE, fill = "#6DCD59FF",
remove.intercepts = FALSE, by = NULL, gsub.pattern = NULL, gsub.replacement = NULL,
xlab = "Effect size", ylab = "Parameter") {
# create objects just to avoid errors with ggplot functions when checking
# package
position_dodge <- org.variable <- value <- significance <- `l-95% CI` <- `u-95% CI` <- theme <- unit <- NULL
if (!is.null(by)) {
levels <- draws[, by]
unique_levels <- unique(levels)
} else unique_levels <- ".A"
# run loop over each level (or just the single 'level')
results_list <- lapply(unique_levels, function(x) {
if (identical(unique_levels, ".A")) {
# keep only betas
draws <- draws[, grep("^b_", names(draws), value = TRUE)]
} else {
draws <- draws[draws[, by, drop = TRUE] == x, grep("^b_", names(draws),
value = TRUE)]
}
# remove intercept betas
if (remove.intercepts)
draws <- draws[, grep("^b_Intercept", names(draws), value = TRUE, invert = TRUE),
drop = FALSE]
# compute model-averaged posteriors of overlapping parameters
coef_table <- posterior::summarise_draws(draws, median, ~quantile(.x, probs = c(0.025,
0.975)), posterior::default_convergence_measures())
names(coef_table)[3:4] <- c("l-95% CI", "u-95% CI")
coef_table$value <- coef_table$median
coef_table$significance <- ifelse(coef_table$`l-95% CI` * coef_table$`u-95% CI` >
0, "sig", "non-sig")
coef_table$significance <- factor(coef_table$significance, levels = c("non-sig",
"sig"))
if (ncol(draws) > 1)
sdraws <- stack(draws[, grep("^b_", names(draws), value = TRUE)], ind = "levels") else sdraws <- data.frame(value = draws[, 1], variable = names(draws))
names(sdraws) <- c("value", "variable")
sdraws$significance <- sapply(sdraws$variable, function(x) coef_table$significance[as.character(coef_table$variable) ==
x][1])
# add level if (!identical(unique_levels, '.A')) { to coef_table
coef_table$.new.column <- x
names(coef_table)[ncol(coef_table)] <- "levels"
# to sdraws
sdraws$.new.column <- x
names(sdraws)[ncol(sdraws)] <- "levels"
# }
output <- list(coef_table = coef_table, sdraws = sdraws)
return(output)
})
# put both results into a single data frame
coef_table <- do.call(rbind, lapply(results_list, "[[", 1))
sdraws <- do.call(rbind, lapply(results_list, "[[", 2))
if (!is.null(gsub.pattern) & !is.null(gsub.replacement)) {
if (length(gsub.pattern) != length(gsub.replacement))
stop2("'gsub.replacement' and 'gsub.pattern' must have the same length")
for (i in 1:length(gsub.pattern)) {
sdraws$variable <- gsub(pattern = gsub.pattern[i], replacement = gsub.replacement[i],
sdraws$variable)
coef_table$variable <- gsub(pattern = gsub.pattern[i], replacement = gsub.replacement[i],
coef_table$variable)
}
}
# duplicate variable column
coef_table <- as.data.frame(coef_table)
coef_table$org.variable <- coef_table$variable
sdraws$org.variable <- sdraws$variable
if (!identical(unique_levels, ".A")) {
coef_table$variable <- paste(coef_table$levels, coef_table$variable, sep = "-")
coef_table <- coef_table[order(coef_table$org.variable), ]
sdraws$variable <- paste(sdraws$levels, sdraws$variable, sep = "-")
}
rownames(coef_table) <- coef_table$variable
fill_df <- data.frame(level = unique_levels, fill = fill)
if (!identical(unique_levels, ".A"))
coef_table$fill_values <- sapply(coef_table$levels, function(x) fill_df$fill[fill_df$level ==
x]) else coef_table$fill_values <- fill[1]
if (highlight) {
coef_table$col_pointrange <- ifelse(coef_table$significance == "non-sig",
"gray", "black")
# coef_table$fill_values <- ifelse(coef_table$significance == 'no-sig',
# grDevices::adjustcolor(col = coef_table$fill_values, alpha.f = 0.5),
# coef_table$fill_values)
sdraws$significance <- sapply(sdraws$variable, function(x) coef_table$significance[coef_table$variable ==
x])
} else {
coef_table$col_pointrange <- rep("black", nrow(coef_table))
sdraws$significance <- 1
}
pd <- position_dodge(width = 0.05)
# creat plots
gg_distributions <- ggplot2::ggplot(data = sdraws, ggplot2::aes(y = org.variable,
x = value, fill = levels, alpha = if (highlight)
significance else NULL)) + ggplot2::geom_vline(xintercept = 0, col = "black", lty = 2) +
ggdist::stat_halfeye(ggplot2::aes(x = value), .width = c(0.95), normalize = "panels",
color = "transparent", position = pd) + ggplot2::scale_alpha_manual(values = c(0.4,
0.8), guide = "none") + ggplot2::scale_fill_manual(values = if (!identical(unique_levels,
".A"))
as.vector(coef_table$fill_values) else fill[1]) + ggplot2::geom_point(data = coef_table, position = pd) + ggplot2::geom_errorbar(data = coef_table,
ggplot2::aes(xmin = `l-95% CI`, xmax = `u-95% CI`), width = 0, position = pd) +
ggplot2::facet_wrap(~org.variable, scales = "free_y", ncol = 1) + ggplot2::theme_classic() +
ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"), plot.margin = ggplot2::margin(0,
0, 0, 0, "pt"), legend.position = "none", strip.background = ggplot2::element_blank(),
strip.text = ggplot2::element_blank()) + ggplot2::labs(x = xlab, y = ylab,
fill = "Effect") + theme(panel.spacing.y = unit(0, "null") # no vertical space between panels
)
coef_table$variable <- coef_table$significance <- coef_table$value <- NULL
names(coef_table) <- c("Estimate", "l-95% CI", "u-95% CI", "Rhat", "Bulk_ESS",
"Tail_ESS")
coef_table <- coef_table[, c("Estimate", "l-95% CI", "u-95% CI", "Rhat", "Bulk_ESS",
"Tail_ESS")]
if (!is.null(name))
cat("\n\n## ", name, "\n\n")
html_coef_table <- brmsish:::html_format_coef_table(coef_table, fill = fill,
highlight = highlight)
print(html_coef_table)
print(gg_distributions)
}1 Read data
1.1 Set working directory
This is the path were the supplementary data is saved and were the output files will be saved. It also creates a directory there in which to save single and averaged models:
Code
# Set working directory
path <- "PUT_WORKING_DIRECTORY_PATH_HERE"
# create directory for single models
dir.create(file.path(path, "single_models"))
# create directory for average models
dir.create(file.path(path, "averaged_models"))The file ‘acoustic_and_climatic_data_by_hour.csv’ contains all the data required to run statistical models down below:
Code
call_rate_hour <- read.csv(file.path(path, "acoustic_and_climatic_data_by_hour.csv"))Describe environmental variables:
Code
agg <- aggregate(cbind(temp, prev_temp, HR, rain, rain_24, rain_48, moonlight) ~
1, call_rate_hour, function(x) round(c(mean(x), sd(x), min(x), max(x)), 3))
agg <- as.data.frame(matrix(unlist(agg), ncol = 4, byrow = TRUE, dimnames = list(c("Temperature",
"Previous temperature", "Relative humidity", "Night rain", "Rain 24 hours", "Rain 48 hours",
"Moonlight"), c("mean", "sd", "min", "max"))))
# print table as kable
kb <- kable(agg, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| mean | sd | min | max | |
|---|---|---|---|---|
| Temperature | 24.229 | 2.126 | 19.472 | 31.763 |
| Previous temperature | 24.097 | 0.988 | 20.005 | 27.111 |
| Relative humidity | 90.056 | 8.464 | 55.158 | 99.940 |
| Night rain | 0.079 | 0.457 | 0.000 | 10.417 |
| Rain 24 hours | 1.441 | 3.161 | 0.000 | 25.620 |
| Rain 48 hours | 1.453 | 3.170 | 0.000 | 25.620 |
| Moonlight | 0.495 | 0.355 | 0.000 | 1.000 |
2 Directed acyclical graphs (DAGs)
The code creates the DAGs used for determining the adjustments sets needed to evaluate causality:
Code
coords <- list(x = c(sc_rain = -0.4, evotranspiration = 0.5, sc_prev_rain = 0.7,
sc_temp = -0.8, sc_HR = 0, n_call = 0, sc_moonlight = 0.3, hour = -0.5), y = c(sc_rain = 0.4,
evotranspiration = 0.3, sc_prev_rain = -0.5, sc_temp = 0, climate = 1, sc_HR = -0.6,
n_call = 0, sc_moonlight = 1, hour = 0.9))
# sc_temp + sc_HR + sc_moonlight + sc_rain + sc_rain_24 + ar(p = 2, time =
# hour_diff, gr = hour sc_temp = temp y meanT = prev_temp
dag_l <- dagify(sc_rain ~ evotranspiration, sc_prev_rain ~ evotranspiration, sc_temp ~
climate, sc_temp ~ sc_rain, sc_HR ~ sc_rain, n_call ~ sc_HR, n_call ~ hour, n_call ~
sc_moonlight, sc_moonlight ~ hour, sc_temp ~ hour, sc_HR ~ sc_temp, sc_HR ~ sc_prev_rain,
sc_HR ~ sc_rain, n_call ~ sc_temp, n_call ~ sc_prev_rain, n_call ~ sc_rain, labels = c(n_call = "Calling\nactivity",
sc_HR = "Relative\nhumidity", sc_rain = "Current\nrain", sc_prev_rain = "Prior\nrain",
sc_moonlight = "Moonlight", hour = "Earth\nrotation", sc_temp = "Tempera-\nture",
evotranspiration = "Evotrans-\npiration", climate = "Climate", latent = c("evotranspiration",
"climate"), outcome = "n_call"), coords = coords)
tidy_dag <- tidy_dagitty(dag_l)
tidy_dag$data$type <- ifelse(is.na(tidy_dag$data$to), "outcome", "predictor")
tidy_dag$data$type[tidy_dag$data$name %in% c("evotranspiration", "climate")] <- "latent"
dat <- tidy_dag$data
shorten_distance <- c(0.07, 0.07)
dat$slope <- (dat$yend - dat$y)/(dat$xend - dat$x)
distance <- sqrt((dat$xend - dat$x)^2 + (dat$yend - dat$y)^2)
proportion <- shorten_distance[1]/distance
dat$xend <- (1 - proportion/2) * dat$xend + (proportion/2 * dat$x)
dat$yend <- (1 - proportion/2) * dat$yend + (proportion/2 * dat$y)
proportion <- shorten_distance[2]/distance
dat$xstart <- (1 - proportion/2) * (dat$x - dat$xend) + dat$xend
dat$ystart <- (1 - proportion/2) * (dat$y - dat$yend) + dat$yend
tidy_dag$data <- dat
basic_dag <- ggplot(tidy_dag, aes(x = x, y = y, xend = xend, yend = yend)) + scale_color_viridis_d(begin = 0.2,
end = 0.8, alpha = 0.5) + geom_dag_text(color = "black", aes(label = label, color = label)) +
labs(color = "Type") + theme_dag() + theme(legend.position = "bottom") + guides(colour = guide_legend(override.aes = list(size = 10))) +
geom_dag_point(aes(color = type), size = 30) + expand_limits(y = c(-0.67, 1.1)) +
geom_dag_edges_fan(edge_color = viridis(10, alpha = 0.4)[2], arrow = grid::arrow(length = grid::unit(10,
"pt"), type = "closed"), aes(x = xstart, y = ystart, xend = xend, yend = yend))
basic_dagCode
dag_24 <- dagify(sc_rain ~ evotranspiration, sc_rain_24 ~ evotranspiration, sc_temp ~
climate, sc_temp ~ sc_rain, sc_HR ~ sc_rain, n_call ~ sc_HR, n_call ~ hour, n_call ~
sc_moonlight, sc_moonlight ~ hour, sc_temp ~ hour, sc_HR ~ sc_temp, sc_HR ~ sc_rain_24,
sc_HR ~ sc_rain, n_call ~ sc_temp, n_call ~ sc_rain_24, n_call ~ sc_rain, labels = c(n_call = "Call rate",
sc_HR = "Relative humidity", sc_rain = "Night Rain", sc_rain_24 = "Previous Rain",
sc_moonlight = "Moonlight", hour = "Earth rotation", sc_temp = "Temperature",
evotranspiration = "Evotranspiration", climate = "Climate", latent = c("evotranspiration",
"climate"), outcome = "n_call"))
dag_48 <- dagify(sc_rain ~ evotranspiration, sc_rain_48 ~ evotranspiration, sc_temp ~
climate, sc_temp ~ sc_rain, sc_HR ~ sc_rain, n_call ~ sc_HR, n_call ~ hour, n_call ~
sc_moonlight, sc_moonlight ~ hour, sc_temp ~ hour, sc_HR ~ sc_temp, sc_HR ~ sc_rain_48,
sc_HR ~ sc_rain, n_call ~ sc_temp, n_call ~ sc_rain_48, n_call ~ sc_rain, labels = c(n_call = "Call rate",
sc_HR = "Relative humidity", sc_rain = "Night Rain", sc_rain_48 = "Previous Rain",
sc_moonlight = "Moonlight", hour = "Earth rotation", sc_temp = "Temperature",
evotranspiration = "Evotranspiration", climate = "Climate", latent = c("evotranspiration",
"climate"), outcome = "n_call"))3 Bayesian regression models
3.1 Scale variables and set model parameters
Prepare data for regression models:
Code
# make hour a factor
call_rate_hour$hour <- factor(call_rate_hour$hour)
# scale and mean-center
call_rate_hour$sc_temp <- scale(call_rate_hour$temp)
call_rate_hour$sc_HR <- scale(call_rate_hour$HR)
call_rate_hour$sc_rain <- scale(call_rate_hour$rain)
call_rate_hour$sc_rain_24 <- scale(call_rate_hour$rain_24)
call_rate_hour$sc_rain_48 <- scale(call_rate_hour$rain_48)
call_rate_hour$sc_moonlight <- scale(call_rate_hour$moonlight)
priors <- c(prior(normal(0, 4), class = "b"))
chains <- 4
iter <- 100003.2 Fit all models
This code fits all the models in the adjustment sets across all predictors:
Code
param_grid <- expand.grid(dag = c("dag_24", "dag_48"), exposure = c("sc_temp", "sc_HR",
"sc_moonlight", "sc_rain", "sc_rain_24", "sc_rain_48"), effect = c("total", "direct"),
stringsAsFactors = FALSE)
param_grid$name <- apply(param_grid, 1, paste, collapse = "-")
# remove wrong dags for previous rain
param_grid <- param_grid[!(param_grid$dag == "dag_24" & param_grid$exposure == "sc_rain_48") &
!(param_grid$dag == "dag_48" & param_grid$exposure == "sc_rain_24"), ]
adjustment_sets_list <- pblapply(seq_len(nrow(param_grid)), cl = 1, function(x) {
forms <- adjustment_set_formulas(dag = if (param_grid$dag[x] == "dag_24")
dag_24 else dag_48, type = if (param_grid$effect[x] == "total")
"all" else "minimal", exposure = param_grid$exposure[x], outcome = "n_call", effect = param_grid$effect[x],
required_variable = "hour", formula_parts = c("n_call | resp_rate(rec_time)",
"+ ar(p = 2, time = hour_diff, gr = hour)"), latent = c("evotranspiration",
"climate"), remove = "hour", plot = FALSE)
return(forms)
})
names(adjustment_sets_list) <- param_grid$name
param_grid$model <- sapply(seq_len(nrow(param_grid)), function(x) {
if (param_grid$effect[x] == "direct")
adjustment_sets_list[[which(names(adjustment_sets_list) == param_grid$name[x])]] else NA
})
param_grid$model <- unlist(param_grid$model)
param_grid <- as.data.frame(param_grid)
param_grid$model[!is.na(param_grid$model)] <- remove_special_chars(param_grid$model[!is.na(param_grid$model)])
param_grid$model <- c("total_effect_temperature_with_rain_24", "total_effect_temperature_with_rain_48",
"total_effect_humidity_with_rain_24", "total_effect_humidity_with_rain_48", "total_effect_moon_with_rain_24",
"total_effect_moon_with_rain_48", "total_effect_rain_with_rain_24", "total_effect_rain_with_rain_48",
"total_effect_previous_rain_24", "total_effect_previous_rain_48", param_grid$model[!is.na(param_grid$model)])
param_grid$exposure.name <- param_grid$exposure
param_grid$exposure.name[grep("temp", param_grid$exposure.name)] <- "Temperature"
param_grid$exposure.name[grep("HR", param_grid$exposure.name)] <- "Relative humidity"
param_grid$exposure.name[grep("moon", param_grid$exposure.name)] <- "Moonlight"
param_grid$exposure.name[grep("rain$", param_grid$exposure.name)] <- "Current rain"
param_grid$exposure.name[grep("rain_24", param_grid$exposure.name)] <- "Previous rain (24h)"
param_grid$exposure.name[grep("rain_48", param_grid$exposure.name)] <- "Previous rain (48h)"
table(param_grid$exposure.name)
write.csv(x = param_grid, file = file.path(path, "direct_and_total_effect_model_data_frame.csv"),
row.names = FALSE)3.3 Combined models need to infer causality
The code below takes all models representing an adjustment set for each predictor and average them into a single model fit:
Code
direct_adjustment_sets_list <- adjustment_sets_list[grep("direct", names(adjustment_sets_list))]
for(i in seq_along(direct_adjustment_sets_list))
pa_comb_mod <-
averaged_model(
formulas = direct_adjustment_sets_list[[i]],
data = call_rate_hour,
suffix = "direct",
model_call = "brm(formula, data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(), prior = priors, backend = 'cmdstanr')",
save.files = TRUE,
path = file.path(path, "averaged_models"),
# name = "temperature_with_rain_24",
cores = 1
)
model_call = "brm(formula, data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(), prior = priors, backend = 'cmdstanr')"
formulas <- unlist(direct_adjustment_sets_list)
mod_path <- file.path(path, "single_models")
fit_list <- pblapply_brmsish_int(X = formulas, cl = 1, function(y) {
# make file name without special characters
mod_name <-
paste0(remove_special_chars(as.character(y)), ".RDS")
if (!file.exists(file.path(mod_path, mod_name))) {
cat("Fitting", y, "\n")
mc <-
gsub(pattern = "formula",
replacement = as.character(y),
x = model_call)
mc <- parse(text = mc)
fit <- eval(mc)
if (save.files)
saveRDS(fit, file = file.path(mod_path, mod_name))
}
})3.4 Results
Code
param_grid <- read.csv(file = file.path(path, "direct_and_total_effect_model_data_frame.csv"))
param_grid$files <- file.path(file.path(path, "averaged_models", paste0(param_grid$model,
".RDS")))
for (i in unique(param_grid$exposure.name)) {
Y <- param_grid[param_grid$exposure.name == i, ]
cat(paste("\n###", i), "\n")
cat("\n#### Direct effects\n")
for (e in which(Y$effect == "direct")) {
if (grepl("24", Y$model[e]))
cat("\n##### 24 hour previous rain model:\n") else cat("\n##### 48 hour previous rain model:\n")
extended_summary(read.file = Y$files[e], highlight = TRUE, remove.intercepts = TRUE,
print.name = FALSE)
cat("\n")
}
cat("\n#### Total effect\n")
for (w in which(Y$effect == "total")) {
if (grepl("24", Y$files[w]))
cat("\n##### 24 hour previous rain model:\n") else cat("\n##### 48 hour previous rain model:\n")
draws <- readRDS(Y$files[w])
draw_extended_summary(draws, highlight = TRUE, remove.intercepts = TRUE)
cat("\n")
cat("\n###### Summary of single models:\n")
# print summary
print(readRDS(gsub("\\.RDS", "_fit_summary.RDS", Y$files[w])))
}
cat("\n")
}3.4.1 Temperature
3.4.1.1 Direct effects
3.4.1.1.1 24 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 13808 | 14519.5 | 1687235322 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_temp | 0.577 | 0.498 | 0.656 | 1 | 13808.0 | 14685.6 |
| b_sc_HR | 0.314 | 0.242 | 0.387 | 1 | 18992.0 | 14519.5 |
| b_sc_rain | 0.037 | 0.002 | 0.074 | 1 | 23964.4 | 15266.3 |
| b_sc_rain_24 | 0.096 | 0.057 | 0.135 | 1 | 23319.6 | 15914.6 |
3.4.1.1.2 48 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 9460.49 | 12432.4 | 293904519 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_temp | 0.565 | 0.487 | 0.643 | 1 | 9460.49 | 12432.4 |
| b_sc_HR | 0.307 | 0.235 | 0.379 | 1 | 11994.69 | 13312.6 |
| b_sc_rain | 0.029 | -0.007 | 0.066 | 1 | 17566.34 | 16050.0 |
| b_sc_rain_48 | 0.006 | -0.032 | 0.044 | 1 | 17289.44 | 16222.0 |
3.4.1.2 Total effect
3.4.1.2.1 24 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_temp | 0.308 | 0.258 | 0.361 | 1.043 | 44.190 | 893.839 |
| b_sc_rain | 0.041 | 0.004 | 0.078 | 1.026 | 152.336 | 830.415 |
3.4.1.2.1.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2542.43 | 5383.76 | 629626678 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2238.34 | 5028.64 | 2096770826 |
| 3 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2320.90 | 4619.08 | 1430954384 |
| 4 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2292.54 | 4473.95 | 1881604582 |
3.4.1.2.2 48 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_temp | 0.305 | 0.258 | 0.353 | 1.005 | 961.292 | 983.283 |
| b_sc_rain | 0.036 | -0.002 | 0.076 | 1 | 1113.522 | 969.564 |
3.4.1.2.2.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2529.79 | 5016.37 | 1471048663 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2588.68 | 5844.63 | 412112948 |
| 3 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2622.63 | 5641.91 | 41093237 |
| 4 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2334.61 | 4741.38 | 366114719 |
3.4.2 Relative humidity
3.4.2.1 Direct effects
3.4.2.1.1 24 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 11039.8 | 13170.1 | 1554924087 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_HR | 0.314 | 0.240 | 0.388 | 1 | 13006.8 | 13679.9 |
| b_sc_rain | 0.037 | 0.002 | 0.075 | 1 | 18235.1 | 15400.0 |
| b_sc_rain_24 | 0.096 | 0.057 | 0.135 | 1 | 16966.8 | 15476.0 |
| b_sc_temp | 0.578 | 0.499 | 0.657 | 1 | 11039.8 | 13170.1 |
3.4.2.1.2 48 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 8448.99 | 11500.1 | 1931567299 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_HR | 0.308 | 0.234 | 0.380 | 1 | 10649.49 | 12433.6 |
| b_sc_rain | 0.029 | -0.008 | 0.067 | 1 | 16664.06 | 15531.5 |
| b_sc_rain_48 | 0.005 | -0.032 | 0.043 | 1 | 15745.61 | 14835.4 |
| b_sc_temp | 0.566 | 0.487 | 0.645 | 1 | 8448.99 | 11500.1 |
3.4.2.2 Total effect
3.4.2.2.1 24 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_HR | 0.309 | 0.240 | 0.383 | 1.006 | 942.064 | 903.519 |
| b_sc_rain | 0.038 | 0.002 | 0.076 | 1.003 | 995.090 | 912.254 |
| b_sc_rain_24 | 0.094 | 0.057 | 0.133 | 1.003 | 993.441 | 772.255 |
| b_sc_temp | 0.575 | 0.496 | 0.655 | 1.004 | 1022.397 | 861.998 |
3.4.2.2.1.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_HR + sc_moonlight + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2556.01 | 5127.35 | 1982612331 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2182.06 | 4328.43 | 1554924087 |
3.4.2.2.2 48 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_HR | 0.303 | 0.230 | 0.371 | 1.006 | 801.527 | 931.946 |
| b_sc_rain | 0.030 | -0.008 | 0.066 | 1.002 | 812.238 | 1071.444 |
| b_sc_rain_48 | 0.003 | -0.034 | 0.044 | 1.005 | 889.913 | 984.434 |
| b_sc_temp | 0.564 | 0.480 | 0.644 | 1.003 | 821.550 | 773.876 |
3.4.2.2.2.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_HR + sc_moonlight + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2326.55 | 4696.55 | 846435536 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2166.28 | 4627.63 | 1931567299 |
3.4.3 Moonlight
3.4.3.1 Direct effects
3.4.3.1.1 48 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 30283.5 | 15929.8 | 1624275037 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_moonlight | -0.129 | -0.193 | -0.064 | 1 | 30283.5 | 15929.8 |
3.4.3.1.2 48 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 30283.5 | 15929.8 | 1624275037 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_moonlight | -0.129 | -0.193 | -0.064 | 1 | 30283.5 | 15929.8 |
3.4.3.2 Total effect
3.4.3.2.1 24 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_moonlight | -0.101 | -0.17 | -0.039 | 1.04 | 40.665 | 750.448 |
3.4.3.2.1.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2463.32 | 4856.49 | 1624275037 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2613.63 | 5435.51 | 1331671454 |
| 3 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2426.99 | 4881.26 | 1714115616 |
| 4 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2184.13 | 4557.11 | 946500998 |
| 5 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2386.65 | 5083.35 | 2088023938 |
| 6 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2571.23 | 4971.80 | 1587538237 |
| 7 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2319.91 | 4626.15 | 1236984196 |
| 8 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2197.61 | 4746.52 | 1987527846 |
| 9 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2353.89 | 5263.44 | 924618251 |
| 10 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2409.47 | 4652.12 | 1063522846 |
| 11 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2407.64 | 4851.64 | 259660598 |
| 12 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2473.17 | 4731.37 | 1625659188 |
| 13 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2596.10 | 5292.48 | 1259975215 |
| 14 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 1998.35 | 4047.25 | 192228769 |
| 15 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2365.76 | 4743.96 | 559459570 |
| 16 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2432.09 | 4364.34 | 1786523665 |
3.4.3.2.2 48 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_moonlight | -0.105 | -0.166 | -0.038 | 1.02 | 111.506 | 659.054 |
3.4.3.2.2.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2534.83 | 5170.40 | 169376128 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2664.91 | 5585.82 | 731739252 |
| 3 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2449.93 | 4958.79 | 230289957 |
| 4 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2283.15 | 4266.04 | 1735286727 |
| 5 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2486.75 | 4561.76 | 1539258099 |
| 6 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2525.61 | 5053.20 | 1234320719 |
| 7 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2481.99 | 4862.29 | 183236870 |
| 8 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2296.24 | 4114.83 | 239469808 |
| 9 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2298.87 | 4495.74 | 1126617762 |
| 10 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2442.38 | 4594.40 | 684931243 |
| 11 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 1 (5e-05%) | 0 | 2416.00 | 4930.19 | 1093486403 |
| 12 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2146.30 | 4395.84 | 1670597608 |
| 13 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2450.19 | 4442.71 | 627408872 |
| 14 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2513.43 | 4849.68 | 1827601939 |
| 15 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2101.51 | 4415.14 | 104020182 |
| 16 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2495.11 | 5464.93 | 598459253 |
| 17 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2120.85 | 4487.97 | 2041361978 |
| 18 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2336.51 | 4816.06 | 1759350038 |
| 19 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 1 (5e-05%) | 0 | 2453.12 | 5040.07 | 927628753 |
| 20 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2360.45 | 5250.70 | 2028572807 |
| 21 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2642.43 | 5077.10 | 293339431 |
| 22 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2342.88 | 4887.79 | 1868346384 |
| 23 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2515.33 | 5176.64 | 552665161 |
| 24 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2253.53 | 4141.73 | 1764642823 |
| 25 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2634.31 | 5119.83 | 505424120 |
| 26 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2427.59 | 5073.57 | 1761124445 |
| 27 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2609.57 | 5134.76 | 1073467634 |
| 28 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2461.43 | 4814.97 | 1419285187 |
| 29 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2596.36 | 5583.40 | 1232772932 |
| 30 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2271.50 | 5053.09 | 278133969 |
| 31 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 1 (5e-05%) | 0 | 2507.89 | 4315.91 | 866725413 |
| 32 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2126.46 | 4549.18 | 1885323046 |
3.4.4 Current rain
3.4.4.1 Direct effects
3.4.4.1.1 24 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain + sc_HR + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 10228.7 | 12456.7 | 1281527469 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain | 0.037 | 0.001 | 0.075 | 1 | 18004.9 | 15909.5 |
| b_sc_HR | 0.313 | 0.240 | 0.386 | 1 | 13025.4 | 13472.6 |
| b_sc_rain_24 | 0.096 | 0.059 | 0.134 | 1 | 17444.9 | 15779.0 |
| b_sc_temp | 0.577 | 0.499 | 0.655 | 1 | 10228.7 | 12456.7 |
3.4.4.1.2 48 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain + sc_HR + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 7371.36 | 11882.6 | 1129340659 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain | 0.029 | -0.007 | 0.066 | 1 | 14353.91 | 14479.2 |
| b_sc_HR | 0.307 | 0.235 | 0.380 | 1 | 8890.14 | 12279.5 |
| b_sc_rain_48 | 0.005 | -0.032 | 0.043 | 1 | 13843.68 | 14097.9 |
| b_sc_temp | 0.566 | 0.487 | 0.644 | 1 | 7371.36 | 11882.6 |
3.4.4.2 Total effect
3.4.4.2.1 24 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain | 0.000 | -0.034 | 0.035 | 1.001 | 942.79 | 942.018 |
| b_sc_rain_24 | 0.076 | 0.039 | 0.113 | 1 | 1056.52 | 982.799 |
3.4.4.2.1.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain + sc_moonlight + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2531.75 | 4835.39 | 1596946440 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2369.16 | 4574.59 | 1687235322 |
3.4.4.2.2 48 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain | -0.005 | -0.042 | 0.030 | 1 | 1024.455 | 819.064 |
| b_sc_rain_48 | 0.011 | -0.027 | 0.049 | 1.009 | 932.089 | 975.233 |
3.4.4.2.2.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain + sc_moonlight + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2647.66 | 5129.91 | 1485832077 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2675.26 | 5099.07 | 167092432 |
3.4.5 Previous rain (24h)
3.4.5.1 Direct effects
3.4.5.1.1 24 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 11374.2 | 13445.2 | 1379661434 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain_24 | 0.096 | 0.058 | 0.135 | 1 | 20345.7 | 15663.3 |
| b_sc_HR | 0.314 | 0.240 | 0.386 | 1 | 14938.4 | 13445.2 |
| b_sc_rain | 0.038 | 0.001 | 0.074 | 1 | 21502.1 | 16005.3 |
| b_sc_temp | 0.578 | 0.498 | 0.656 | 1 | 11374.2 | 13493.9 |
3.4.5.2 Total effect
3.4.5.2.1 24 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain_24 | 0.088 | 0.051 | 0.127 | 1.01 | 883.631 | 733.445 |
| b_sc_rain | 0.043 | 0.005 | 0.081 | 1.001 | 908.203 | 446.195 |
3.4.5.2.1.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2692.30 | 5204.02 | 423216939 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2556.12 | 4635.90 | 610868624 |
| 3 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 1 (5e-05%) | 0 | 2559.09 | 4449.50 | 1988187279 |
| 4 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2077.52 | 4134.53 | 640249946 |
3.4.6 Previous rain (48h)
3.4.6.1 Direct effects
3.4.6.1.1 48 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain_48 + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 11112 | 13056.2 | 1238572765 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain_48 | 0.006 | -0.031 | 0.043 | 1 | 17211.0 | 15210.0 |
| b_sc_HR | 0.308 | 0.236 | 0.380 | 1 | 14154.7 | 13956.7 |
| b_sc_rain | 0.029 | -0.007 | 0.066 | 1 | 19170.2 | 15964.2 |
| b_sc_temp | 0.566 | 0.486 | 0.645 | 1 | 11112.0 | 13056.2 |
3.4.6.2 Total effect
3.4.6.2.1 48 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain_48 | 0.005 | -0.033 | 0.043 | 1.014 | 711.511 | 873.701 |
| b_sc_rain | 0.036 | 0.001 | 0.073 | 1 | 857.372 | 792.928 |
3.4.6.2.1.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2447.89 | 5823.08 | 1081516016 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 2 (1e-04%) | 0 | 2413.93 | 4851.86 | 217409394 |
| 3 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2506.36 | 5173.40 | 724206553 |
| 4 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2641.41 | 5493.32 | 509251240 |
| 5 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_48 + sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2386.41 | 4640.41 | 518390312 |
| 6 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_48 + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2575.04 | 5332.87 | 76427694 |
3.5 Combined results with causal inference estimates
Takes the posterior probability distributions from the right causal models
Code
param_grid <- read.csv(file = file.path(path, "direct_and_total_effect_model_data_frame.csv"))
param_grid <- param_grid[param_grid$effect == "direct", ]
param_grid$file <- paste0(remove_special_chars(param_grid$model), ".RDS")
rdss_24 <- list.files(file.path(path, "averaged_models"), pattern = "24.RDS", full.names = TRUE)
combined_draws_list <- lapply(rdss_24, function(x) {
total_draws <- readRDS(x)
exp <- attributes((attributes(total_draws)$averaged_fit_formulas))$exposure.formula
exp <- gsub("n_call | resp_rate(rec_time) ~ ", "", exp, fixed = TRUE)
exposure <- exp <- gsub(" + ar(p = 2, time = hour_diff, gr = hour)", "", exp,
fixed = TRUE)
exp <- paste0("b_", exp)
total_draws <- total_draws[, colnames(total_draws) == exp, drop = FALSE]
names(total_draws) <- exp
direct_fit_file <- param_grid$file[param_grid$exposure == exposure]
direct_fit_file <- direct_fit_file[!duplicated(direct_fit_file)]
if (length(direct_fit_file) > 1)
direct_fit_file <- grep("24", direct_fit_file, value = TRUE)
direct_fit <- readRDS(file = file.path(path, "averaged_models", direct_fit_file))
direct_draws <- posterior::merge_chains(as_draws(direct_fit, variable = exp))
direct_draws <- as.data.frame(thin_draws(direct_draws, thin = length(direct_draws[[1]][[1]])/(nrow(total_draws)))[[1]])
direct_draws$effect <- "direct"
total_draws$effect <- "total"
draws <- rbind(direct_draws, total_draws)
return(draws)
})
combined_draws <- do.call(cbind, combined_draws_list)
combined_draws <- combined_draws[, c(which(sapply(combined_draws, is.numeric)), ncol(combined_draws))]
combined_draws[, -ncol(combined_draws)] <- to_change_percentage(combined_draws[,
-ncol(combined_draws)])
# combined_draws <- as.data.frame(combined_draws)
combined_draws$effect <- ifelse(combined_draws$effect == "direct", "Direct", "Total")
saveRDS(combined_draws, file.path(path, "combined_draws_for_total_and_direct_effects_24h_previous_rain.RDS"))Code
combined_draws <- readRDS(file.path(path, "combined_draws_for_total_and_direct_effects_24h_previous_rain.RDS"))
fill_colors <- viridis::mako(10)[c(8, 4)]
gg_dists <- draw_extended_summary(draws = combined_draws, highlight = TRUE, remove.intercepts = TRUE,
fill = adjustcolor(fill_colors, alpha.f = 0.4), by = "effect", gsub.pattern = c("b_sc_HR",
"b_sc_rain$", "b_sc_rain_24", "b_sc_temp", "b_sc_moonlight"), gsub.replacement = c("Relative\nhumidity",
"Current\nrain", "Prior\nrain", "Temperature", "Moonlight"), ylab = "Variable",
xlab = "Effect size (% of change)")| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| Direct-Current rain | 3.652 | -0.002 | 7.942 | 1.009 | 999.221 | 949.102 |
| Total-Current rain | 0.000 | -3.373 | 3.596 | 1.001 | 942.790 | 942.018 |
| Direct-Moonlight | -12.018 | -17.074 | -6.217 | 1.001 | 972.135 | 944.341 |
| Total-Moonlight | -9.600 | -15.624 | -3.860 | 1.04 | 40.665 | 750.448 |
| Direct-Prior rain | 10.060 | 5.793 | 14.735 | 1 | 916.714 | 983.283 |
| Total-Prior rain | 9.189 | 5.189 | 13.581 | 1.01 | 883.631 | 733.445 |
| Direct-Relative humidity | 37.000 | 27.013 | 47.025 | 1.001 | 902.884 | 848.891 |
| Total-Relative humidity | 36.187 | 27.104 | 46.613 | 1.006 | 942.064 | 903.519 |
| Direct-Temperature | 77.591 | 65.145 | 92.509 | 1.002 | 955.952 | 875.689 |
| Total-Temperature | 36.112 | 29.497 | 43.428 | 1.043 | 44.190 | 893.839 |
Code
gg_dists + ggplot2::scale_fill_manual(values = fill_colors) + ggplot2::theme(axis.ticks.length = ggplot2::unit(0,
"pt"), plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "inside",
legend.position.inside = c(0.7, 0.7))Code
param_grid <- read.csv(file = file.path(path, "direct_and_total_effect_model_data_frame.csv"))
param_grid <- param_grid[param_grid$effect == "direct", ]
param_grid$file <- paste0(remove_special_chars(param_grid$model), ".RDS")
rdss_48 <- list.files(file.path(path, "averaged_models"), pattern = "48.RDS", full.names = TRUE)
combined_draws_list <- lapply(rdss_48, function(x) {
total_draws <- readRDS(x)
exp <- attributes((attributes(total_draws)$averaged_fit_formulas))$exposure.formula
exp <- gsub("n_call | resp_rate(rec_time) ~ ", "", exp, fixed = TRUE)
exposure <- exp <- gsub(" + ar(p = 2, time = hour_diff, gr = hour)", "", exp,
fixed = TRUE)
exp <- paste0("b_", exp)
total_draws <- total_draws[, colnames(total_draws) == exp, drop = FALSE]
names(total_draws) <- exp
direct_fit_file <- param_grid$file[param_grid$exposure == exposure]
direct_fit_file <- direct_fit_file[!duplicated(direct_fit_file)]
if (length(direct_fit_file) > 1)
direct_fit_file <- grep("48", direct_fit_file, value = TRUE)
direct_fit <- readRDS(file = file.path(path, "averaged_models", direct_fit_file))
direct_draws <- posterior::merge_chains(as_draws(direct_fit, variable = exp))
direct_draws <- as.data.frame(thin_draws(direct_draws, thin = length(direct_draws[[1]][[1]])/(nrow(total_draws)))[[1]])
direct_draws$effect <- "direct"
total_draws$effect <- "total"
draws <- rbind(direct_draws, total_draws)
return(draws)
})
combined_draws <- do.call(cbind, combined_draws_list)
combined_draws <- combined_draws[, c(which(sapply(combined_draws, is.numeric)), ncol(combined_draws))]
combined_draws[, -ncol(combined_draws)] <- to_change_percentage(combined_draws[,
-ncol(combined_draws)])
# combined_draws <- as.data.frame(combined_draws)
combined_draws$effect <- ifelse(combined_draws$effect == "direct", "Direct", "Total")
saveRDS(combined_draws, file.path(path, "combined_draws_for_total_and_direct_effects_48h_previous_rain.RDS"))Code
combined_draws <- readRDS(file.path(path, "combined_draws_for_total_and_direct_effects_48h_previous_rain.RDS"))
gg_dists <- draw_extended_summary(draws = combined_draws, highlight = TRUE, remove.intercepts = TRUE,
fill = adjustcolor(fill_colors, alpha.f = 0.4), by = "effect", gsub.pattern = c("b_sc_HR",
"b_sc_rain$", "b_sc_rain_48", "b_sc_temp", "b_sc_moonlight"), gsub.replacement = c("Relative\nhumidity",
"Current\nrain", "Prior\nrain", "Temperature", "Moonlight"), ylab = "Variable",
xlab = "Effect size (% of change)")| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| Direct-Current rain | 2.940 | -0.477 | 6.832 | 1 | 933.292 | 847.036 |
| Total-Current rain | -0.544 | -4.093 | 3.071 | 1 | 1024.455 | 819.064 |
| Direct-Moonlight | -12.018 | -17.074 | -6.217 | 1.001 | 972.135 | 944.341 |
| Total-Moonlight | -9.925 | -15.276 | -3.765 | 1.02 | 111.506 | 659.054 |
| Direct-Prior rain | 0.623 | -2.979 | 4.124 | 1.003 | 857.680 | 1009.747 |
| Total-Prior rain | 0.540 | -3.268 | 4.427 | 1.014 | 711.511 | 873.701 |
| Direct-Relative humidity | 36.477 | 26.300 | 46.337 | 0.999 | 1043.482 | 878.337 |
| Total-Relative humidity | 35.405 | 25.805 | 44.967 | 1.006 | 801.527 | 931.946 |
| Direct-Temperature | 76.199 | 63.368 | 89.305 | 1.004 | 1009.050 | 835.503 |
| Total-Temperature | 35.603 | 29.387 | 42.312 | 1.005 | 961.292 | 983.283 |
Code
gg_dists + ggplot2::scale_fill_manual(values = fill_colors) + ggplot2::theme(axis.ticks.length = ggplot2::unit(0,
"pt"), plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "inside",
legend.position.inside = c(0.7, 0.7))Takeaways
- Variation in call activity strongly linked to environmental variation
Session information
R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: America/Costa_Rica
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] pbapply_1.7-2 tidybayes_3.0.6 ggdag_0.2.12 dagitty_0.3-4
[5] ohun_1.0.2 warbleR_1.1.30 NatureSounds_1.0.4 seewave_2.2.3
[9] tuneR_1.4.6 brmsish_1.0.0 lunar_0.2-1 ggplot2_3.5.1
[13] ggdist_3.3.2 knitr_1.47 kableExtra_1.4.0 HDInterval_0.2.4
[17] readxl_1.4.3 posterior_1.5.0 cowplot_1.1.3 brms_2.21.0
[21] Rcpp_1.0.12 viridis_0.6.5 viridisLite_0.4.2 remotes_2.5.0
[25] sketchy_1.0.3
loaded via a namespace (and not attached):
[1] tensorA_0.36.2.1 rstudioapi_0.16.0 jsonlite_1.8.8
[4] magrittr_2.0.3 TH.data_1.1-2 estimability_1.4.1
[7] farver_2.1.2 rmarkdown_2.27 vctrs_0.6.5
[10] memoise_2.0.1 RCurl_1.98-1.14 htmltools_0.5.8.1
[13] distributional_0.4.0 curl_5.2.0 signal_1.8-0
[16] cellranger_1.1.0 StanHeaders_2.32.9 KernSmooth_2.23-20
[19] htmlwidgets_1.6.4 plyr_1.8.9 sandwich_3.1-0
[22] testthat_3.2.1.1 cachem_1.1.0 emmeans_1.9.0
[25] zoo_1.8-12 igraph_2.0.3 lifecycle_1.0.4
[28] pkgconfig_2.0.3 Matrix_1.6-5 R6_2.5.1
[31] fastmap_1.2.0 digest_0.6.36 colorspace_2.1-0
[34] labeling_0.4.3 fansi_1.0.6 polyclip_1.10-6
[37] abind_1.4-5 compiler_4.3.2 proxy_0.4-27
[40] withr_3.0.0 backports_1.5.0 inline_0.3.19
[43] DBI_1.2.2 QuickJSR_1.2.2 pkgbuild_1.4.4
[46] highr_0.11 ggforce_0.4.2 MASS_7.3-55
[49] rjson_0.2.21 classInt_0.4-10 loo_2.7.0
[52] tools_4.3.2 units_0.8-5 ape_5.8
[55] fftw_1.0-8 glue_1.7.0 nlme_3.1-155
[58] grid_4.3.2 sf_1.0-15 checkmate_2.3.1
[61] reshape2_1.4.4 generics_0.1.3 diffobj_0.3.5
[64] gtable_0.3.5 class_7.3-20 tidyr_1.3.1
[67] tidygraph_1.3.1 xml2_1.3.6 utf8_1.2.4
[70] ggrepel_0.9.5 pillar_1.9.0 stringr_1.5.1
[73] splines_4.3.2 tweenr_2.0.2 dplyr_1.1.4
[76] lattice_0.20-45 survival_3.2-13 tidyselect_1.2.1
[79] arrayhelpers_1.1-0 gridExtra_2.3 V8_4.4.2
[82] svglite_2.1.3 stats4_4.3.2 xfun_0.45
[85] graphlayouts_1.1.1 bridgesampling_1.1-2 brio_1.1.4
[88] matrixStats_1.3.0 rstan_2.32.6 stringi_1.8.4
[91] yaml_2.3.8 boot_1.3-28 xaringanExtra_0.7.0
[94] evaluate_0.24.0 codetools_0.2-18 dtw_1.23-1
[97] ggraph_2.2.1 tibble_3.2.1 cli_3.6.3
[100] RcppParallel_5.1.7 xtable_1.8-4 systemfonts_1.1.0
[103] munsell_0.5.1 coda_0.19-4.1 svUnit_1.0.6
[106] parallel_4.3.2 rstantools_2.4.0 bayesplot_1.11.1
[109] Brobdingnag_1.2-9 bitops_1.0-7 mvtnorm_1.2-5
[112] scales_1.3.0 e1071_1.7-14 packrat_0.9.2
[115] purrr_1.0.2 crayon_1.5.2 rlang_1.1.4
[118] multcomp_1.4-25 formatR_1.14